#include <stdio.h>
#include <stdlib.h>

//Numero de thread que irão rodar em cada bloco CUDA. Deve ser uma potência de 2 para que a operação de redução
//funcione corretamente
#define THREADS_BLOCK 256

//Kernel CUDA que calcula cada iteração de resultados. É chamado várias vezes e retorna Xk+1 dado um Xk e as 
//matrizes A* e B*
__global__ void jacobiKernel (float* Ae, float* Be, float* Xk, float* Xk1, int jOrder)
{
	//Vetor utilizado para cálculo do somatório. É visível para todas as threads de um bloco. Sofrerá redução no final
	//da iteração.
	__shared__ float cache[THREADS_BLOCK];
	
	float temp = 0.0;
	int tid = threadIdx.x;
	int cacheIndex = tid;
	int bId = blockIdx.x;

	//Todas as threads fazem blocos de calculos e avançam até calcularem todo a parte de somatório que é de sua 
	//responsabilidade
	while (tid < jOrder)
	{
		temp += Ae[bId * jOrder + tid] * Xk[tid];
		tid += THREADS_BLOCK;
	}	

	//Cada thread armazena sua parte dos calculos para que sejam reduzidas.
	cache[cacheIndex] = temp; 

	//Aguarda para que todas as threads tenham terminado de calcular suas parciais
	__syncthreads();

	//Faz a redução do somatório. Devido à divisão por 2 da linha abaixo, apenas vetores de tamanho em potências de
	//2 são corretamente reduzidos. Funcionamento da redução pode ser exemplificado como se "dobrassemos" o vetor ao
	//meio, somando cada posição correspondente e repetindo várias vezes, até obtermos apenas um valor.
	int i = THREADS_BLOCK/2;
	while (i != 0) 
	{
		if (cacheIndex < i)
		{
			cache[cacheIndex] += cache[cacheIndex + i];
		}
		//Espera que todas as threads somem as posições correspondentes para que possamos "dobrar" novamente o vetor.
		__syncthreads();
		i /= 2;
	}

	//No final, o valor reduzido do vetor se encontra na posição 0 e então calculamos cada elemento do vetor Xk+1
    if (cacheIndex == 0)
    {
		Xk1[bId] = Be[bId] - cache[0];
	}
}

int main(int argc, char* argv[]) 
{
	//Variáveis do algoritmo. 
	float *A, *Ae, *B, *Be;
	float *Xk, *Xk1, *Diff; 
	float jError;
	float result = 0.0;
	int jOrder, jRowTest, jIteMax, i, j;

	//Ponteiros das variáveis que serão enviadas/lidas para/da GPU
	float *dev_Ae, *dev_Be, *dev_Xk, *dev_Xk1;
	FILE* fpin;

	//Leitura do arquivo de entrada, alocação das matrizes e processamento das matrizes A* e B*
	fpin = fopen(argv[1],"r");

	if (fpin == NULL) 
	{
    	printf("There was an error reading the file!\n");
    	exit(1);
	}	

	fscanf(fpin, "%d", &jOrder);

	B = (float *) malloc(jOrder * sizeof(float));
	Be = (float *) malloc(jOrder * sizeof(float));
	A = (float *) malloc(jOrder * jOrder * sizeof(float));
	Ae = (float *) malloc(jOrder * jOrder * sizeof(float));
	Xk = (float *) malloc(jOrder * sizeof(float));
	Xk1 = (float *) malloc(jOrder * sizeof(float));
	Diff = (float *) malloc(jOrder * sizeof(float)); 

	fscanf(fpin, "%d", &jRowTest);
	fscanf(fpin, "%f", &jError);
	fscanf(fpin, "%d", &jIteMax);  

	printf("Input file was read.\nMatrix Order: %d\nRow Test: %d\nError: %f\nMaximum Number of Iterations: %d\n", jOrder, jRowTest, jError, jIteMax);

	for (i = 0; i < jOrder; i++) 
	{
    	for (j = 0; j < jOrder; j++)
    	{
    		fscanf(fpin, "%f", &A[((i*jOrder) + j)]);
    	}
	}

	for (i = 0; i < jOrder; i++) 
	{
    	fscanf(fpin, "%f", &B[i]);
	}  

	for (i = 0; i < jOrder; i++) 
	{
    	float diagonal = A[((i*jOrder)+i)];
    	Be[i] = B[i]/diagonal;
		Xk[i] = Be[i];
    	for (j = 0; j < jOrder; j++)
    	{
    		if (j == i)
    		{
        		Ae[((i*jOrder) + j)] = 0.0;
    		}
    		else
    		{
        		Ae[((i*jOrder) + j)] = A[((i*jOrder) + j)]/diagonal;
    		}
    	}
	}

	//Aloca e copia as matrizes que serão usadas na GPU
	cudaMalloc((void**)&dev_Ae, jOrder * jOrder * sizeof(float));
	cudaMalloc((void**)&dev_Be, jOrder * sizeof(float));
	cudaMalloc((void**)&dev_Xk, jOrder * sizeof(float));
	cudaMalloc((void**)&dev_Xk1, jOrder * sizeof(float));

	//Copia os valores das matrizes A* e B* para a GPU para que sejam feitos os cálculos
	cudaMemcpy(dev_Ae, Ae, jOrder * jOrder * sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(dev_Be, Be, jOrder * sizeof(float), cudaMemcpyHostToDevice);

	//Iterações do algoritmo. Utiliza o Kernel CUDA para calcular Xk+1, calcula o erro e testa se a execução foi
	//satisfatória
	for (i = 0 ; i < jIteMax ; i++)
	{
		//Copia o atual Xk para a GPU
		cudaMemcpy(dev_Xk, Xk, jOrder * sizeof(float), cudaMemcpyHostToDevice);
		//Chama o Kernel para o cálculo de Xk+1 
		jacobiKernel<<<jOrder,THREADS_BLOCK>>>(dev_Ae, dev_Be, dev_Xk, dev_Xk1, jOrder);
		//Copia Xk+1 da GPU para a memória principal
		cudaMemcpy(Xk1, dev_Xk1, jOrder * sizeof(float), cudaMemcpyDeviceToHost);
		float maxDiff = 0.0;
		float maxXk1 = 0.0;
		float absXk1[5000];
		//Calculo do erro
		for (j = 0 ; j < jOrder ; j++)
		{
			Diff[j] = Xk1[j] - Xk[j]; //Calcula Diff[j]
			Diff[j] = Diff[j] < 0 ? -Diff[j] : Diff[j]; //Valor absoluto de Diff[j]
			absXk1[j] = Xk1[j];
			absXk1[j] = absXk1[j] < 0 ? -absXk1[j] : absXk1[j];		
			maxDiff = Diff[j] > maxDiff ? Diff[j] : maxDiff; //Maior valor do vetor Diff
			maxXk1 = absXk1[j] > maxXk1 ? absXk1[j] : maxXk1; //Maior valor absoluto de Xk+1
			Xk[j] = Xk1[j];
		}
		float Mr = maxDiff / maxXk1; //Calculo final do erro Mr
		//Caso o erro tenha sido satisfeito, encerra execução 
		if (Mr <= jError)
		{
			printf("Criterio de parada atingido:  Mr: %f <= %f\n", Mr, jError);
			break;
		}
	}	
	//Executa o RowTest
	for (j = 0 ; j < jOrder ; j++)
	{
		result += A[((jRowTest * jOrder) + j)] * Xk[j];
	}

	//Mostra os resultados
	printf("Iteracoes: %d\n", i+1);
	printf("RowTest: %d => [%f] =? %f\n",jRowTest, result, B[jRowTest]);

	//Desaloca toda a memória alocada em GPU e memória principal.
	cudaFree(dev_Ae);
	cudaFree(dev_Be);
	cudaFree(dev_Xk1);
	cudaFree(dev_Xk);
	free(A);
	free(B);
	free(Ae);
	free(Be);
	free(Xk);
	free(Xk1);
	fclose(fpin);

	return 0;
}
